SoSe2021
Folienübersicht
Kontinuierliche und kategoriale X: Kovarianzanalyse
Kovarianzanalyse oder ANalysis of COVAriance (ANCOVA)
Ein Hybrid zwischen Regression und ANOVA
- Einige betrachten die ANCOVA als ein “ANOVA”-Modell mit einer kontinuierlichen X Variable (Kovariate) → Fokus liegt auf den Faktorstufen, welche für die Kovariate angepasst wurden (z. B. Quinn & Keough, 2002).
- Andere betrachten es als ein “Regressions”-Modell mit einem kategorialen Prädiktor → Fokus liegt auf der Kovariate, welche für den Faktor angepasst wurde (z. B. Crawley 2005,2007).
- Das typische “Maximal”-Modell beinhaltet die Schätzung eines Achsenabschnitts und einer Steigung (Regressionsteil) für jede Stufe der kategorialen Variable(n) (ANOVA-Teil).
Hybrid zwischen Regression und ANOVA
Gewichtsbeispiel 1
Aus der Regressionsperspektive
\[Weight_{ij}=\alpha_i+\beta_i*Age_j + \epsilon_{ij}\]
Gleichung nach Geschlecht getrennt
\[Weight_{female,j}=\alpha_{female}+\beta_{female}*Age_j + \epsilon_{female,j}\] \[Weight_{male,j}=\alpha_{male}+\beta_{male}*Age_j + \epsilon_{male,j}\]
Hybrid zwischen Regression und ANOVA
Gewichtsbeispiel 2
Mögliche Modelltypen
Gewichtsbeispiel 1
Mögliche Modelltypen
Gewichtsbeispiel 2
Mögliche Modelltypen
Gewichtsbeispiel 3
Signifikanztests
Achsenabschnitte \(\alpha_i\)
- Nullhypothese \(H_0: \alpha_1 = \alpha_2 = ... \alpha_i\)
- Bei 2 Achsenabschnitten: t-Test
- Bei >2 Achsenabschnitten: ANOVA (F-Test)
Steigungen \(\beta_i\)
- Nullhypothese \(H_0: \beta_1 = \beta_2 = ... \beta_i\)
- Bei 2 Steigungen: t-Test
- Bei >2 Steigungen: ANOVA (F-Test)
Beispiel aus DS1: Königslachse
ChinookArg
Logarithmierte Gewichtsunterschiede und Längen-Gewichtsbeziehung
Längen und Gewichte für Königslachse von drei Orten in Argentinien (Daten sind im FSA Paket enthalten):
data("ChinookArg", package = "FSA")
str(ChinookArg)
'data.frame': 112 obs. of 3 variables: $ tl : num 120 115 111 110 110 ... $ w : num 17.9 17.2 16.8 15.8 14.3 13.8 12.8 11.7 12.8 14.8 ... $ loc: Factor w/ 3 levels "Argentina","Petrohue",..: 1 1 1 1 1 1 1 1 1 1 ...
Da für die meisten Arten die Beziehung zwischen dem Gewicht eines Tieres und seiner Länge einer 2-Parameter Powerfunktion (\(W = aL^{b}\)) folgt, wollen wir die Variablen log-transformieren:
ChinookArg$tl_log <- log(ChinookArg$tl) ChinookArg$w_log <- log(ChinookArg$w)
Königslachse
Datenexploration
Wie ist die Längen-Gewichtsbeziehung und gibt es Ortsunterschiede?
Königslachse
Maximalmodell
full_mod <- lm(w_log ~ loc + tl_log + loc:tl_log, ChinookArg) # oder: ~ loc * tl_log summary(full_mod)
Call:
lm(formula = w_log ~ loc + tl_log + loc:tl_log, data = ChinookArg)
Residuals:
Min 1Q Median 3Q Max
-0.58273 -0.18471 -0.00186 0.13088 1.63620
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.6750 1.5904 -4.197 5.64e-05 ***
locPetrohue -2.3957 3.1494 -0.761 0.449
locPuyehue -2.0696 1.6868 -1.227 0.223
tl_log 1.9836 0.3530 5.619 1.56e-07 ***
locPetrohue:tl_log 0.4795 0.6928 0.692 0.490
locPuyehue:tl_log 0.3624 0.3793 0.955 0.342
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3201 on 106 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8924
F-statistic: 185 on 5 and 106 DF, p-value: < 2.2e-16
Kleiner Exkurs:
Dummy-Variablen in R
Wenn erklärende Variablen kategorial sind
Falls x kategorial ist, macht es nicht viel Sinn das Model \(y = a + bx\) zu benutzen, da \(x\) nicht mit \(b\) multipliziert werden kann. Aber R hat eine Hilfskonstruktion:
- Die einzelnen Stufen von kategoriale Variablen werden in sog. Dummy-Variablen jeweils konvertiert.
- Jede Dummy-Variable ist als 0 (FALSE) oder 1 (TRUE) kodiert.
- Anzahl der Dummy-Variablen = Anzahl der Gruppen minus 1
- Alle linearen Modelle passen kategoriale erklärende Variablen mit Hilfe von Dummy-Variablen an.
Demonstration mit model.matrix()
Die Daten
df
y length 1 6.7 M 2 16.6 M 3 12.8 M 4 16.1 M 5 4.2 L 6 1.3 S 7 18.6 L 8 5.0 S 9 8.4 S 10 2.0 S
Die Dummy-Variablen
modelr::model_matrix(df, y ~ length)
# A tibble: 10 x 3
`(Intercept)` lengthM lengthS
<dbl> <dbl> <dbl>
1 1 1 0
2 1 1 0
3 1 1 0
4 1 1 0
5 1 0 0
6 1 0 1
7 1 0 0
8 1 0 1
9 1 0 1
10 1 0 1Wo ist die Längen-Klasse L geblieben?
Demonstration mit model.matrix()
Die Daten
df
y length 1 6.7 M 2 16.6 M 3 12.8 M 4 16.1 M 5 4.2 L 6 1.3 S 7 18.6 L 8 5.0 S 9 8.4 S 10 2.0 S
Die Dummy-Variablen
modelr::model_matrix(df, y ~ length)
# A tibble: 10 x 3
`(Intercept)` lengthM lengthS
<dbl> <dbl> <dbl>
1 1 1 0
2 1 1 0
3 1 1 0
4 1 1 0
5 1 0 0
6 1 0 1
7 1 0 0
8 1 0 1
9 1 0 1
10 1 0 1L ist durch den Achsenschnittpunkt (erste Spalte) dargestellt!
Dummy-Variablen am Beispiel der Königslachse
model.matrix(w_log ~ loc, ChinookArg) %>%
as.data.frame(.) %>%
# Hinzufügen der Original loc Variable
mutate(sampl_loc = ChinookArg$loc) %>%
head()
(Intercept) locPetrohue locPuyehue sampl_loc 1 1 0 0 Argentina 2 1 0 0 Argentina 3 1 0 0 Argentina 4 1 0 0 Argentina 5 1 0 0 Argentina 6 1 0 0 Argentina
Wir können sehen, dass die ersten 5 Lachse in ‘Argentinien’ gefangen wurden und dass die Beobachtungen mit einer 1 für die ‘intercept’ Variable (erste Faktorstufe → hier Argentinien) kodiert und 0 für die Variable locPetrohue und locPuyehue wurden.
Interpretation der Koeffizienten
1
Wie würdest Du diese Koeffizienten interpretieren (wenn also nur der Fangort die erklärende Variable ist)?
chin_mod <- lm(w_log ~ loc, ChinookArg) coef(chin_mod)
(Intercept) locPetrohue locPuyehue 2.25605528 -0.09844412 -1.52993775
Wie lautet die Gleichung für den Fangort?
Interpretation der Koeffizienten
2
Setze die berechneten Koeffizienten in die Regressionsgleichung mit Dummy-Variablen ein:
Zurück zum Maximalmodell mit Interaktion
1
summary(full_mod)
Call:
lm(formula = w_log ~ loc + tl_log + loc:tl_log, data = ChinookArg)
Residuals:
Min 1Q Median 3Q Max
-0.58273 -0.18471 -0.00186 0.13088 1.63620
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.6750 1.5904 -4.197 5.64e-05 ***
locPetrohue -2.3957 3.1494 -0.761 0.449
locPuyehue -2.0696 1.6868 -1.227 0.223
tl_log 1.9836 0.3530 5.619 1.56e-07 ***
locPetrohue:tl_log 0.4795 0.6928 0.692 0.490
locPuyehue:tl_log 0.3624 0.3793 0.955 0.342
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3201 on 106 degrees of freedom
Multiple R-squared: 0.8972, Adjusted R-squared: 0.8924
F-statistic: 185 on 5 and 106 DF, p-value: < 2.2e-16Zurück zum Maximalmodell mit Interaktion
2
Grundgleichung
\[w.log_{ij} = \alpha_{i} + \beta_{i}*tl.log_{j} + \epsilon_{ij}\]
wobei i = Fangort, j = jeder einzelne Lachs
Ortspezifische Gleichungen
\[w.log_{A,j} = \alpha_{A} + \beta_{A}*tl.log_{j} + \epsilon_{A,j}\] \[w.log_{Pe,j} = \alpha_{Pe} + \beta_{Pe}*tl.log_{j} + \epsilon_{Pe,j}\] \[w.log_{Pu,j} = \alpha_{Pu} + \beta_{Pu}*tl.log_{j} + \epsilon_{Pu,j}\]
Berechnung der mittleren Gewichte pro Fangort
Gleichung mit Dummy Variablen
\[w.log = a_A + a_{Pe}*locPetrohue + a_{Pu}*locPuyehue +\\ b_A*tl.log+ b_{Pe}*locPetrohue*tl.log + b_{Pu}*locPuyehue*tl.log\]
Einsetzen der Koeffizienten aus R in die Gleichungen
coef(full_mod)
(Intercept) locPetrohue locPuyehue tl_log
-6.6749719 -2.3957331 -2.0695530 1.9836016
locPetrohue:tl_log locPuyehue:tl_log
0.4794543 0.3624014 Argentina: \(w.log_{A} = -6.67+1.98*tl_log\)
Petrohuye: \(w.log_{Pe}= (-6.67+-2.40)+(1.98+0.48*1)*tl.log = -9.07+2.46*tl.log\)
Puyehuye: \(w.log_{Pu} = (-6.67+-2.07)+(1.98+0.36*1)*tl.log = -8.74+2.34*tl.log\)
Systematische Modellauswahl
1
Anstelle der ‘all subsets selection’ (siehe Bsp. airquality), wollen wir hier etwas systematischer vorgehen und die sog. ‘backward selection’ vornehmen.
‘backward selection’
- Starte mit dem Maximalmodell, welches alle X Variablen und mögliche Interaktionen enthält, und überprüfe die Modellannahmen anhand dieses vollen Modells.
- Erstelle nun mehrere reduzierte Modelle, bei der jeweils 1 Interaktionsterm oder ein Hauptterm (der nicht in einem Interaktionsterm vorkommt) entfernt wurde.
- Vergleiche das volle mit den reduzierten Modellen mittels z.B. AIC und wähle das Modell, welches am meisten erklären kann bzw. den niedrigsten AIC hat: dieses stellt nun das neue ‘volle’ Modell dar.
- Erstelle erneut reduzierte Modelle und vergleiche mittels AIC.
- Wiederhole diesen Prozess bis das ‘volle’ Modell den niedrigsten AIC hat.
Systematische Modellauswahl
2
Alternativen
- ‘forward selection’: hier wird genau umgekehrt verfahren, es wird mit dem Nullmodell gestartet und nacheinander eine Variable hinzugefügt.
- ‘stepwise forward/backward selection’:
- Kombiniert beide Ansätze.
R Funktionen
step()→ führt eine vollständige Selektion durch, bis das Modell mit dem niedrigsten AIC gefunden wird.- Mit dem Argument
methodkann eingestellt werden, ob die Richtung ‘forward’, ‘backward’ oder beides (‘both’) sein soll. - Hinweis: Die Funktion wendet nicht die Differenzregel (von mind. 2 AIC Werten) für komplexere Modelle an, d.h. es muss hier evt. manuell weiter selektiert werden!
- Mit dem Argument
drop1→ erstellt reduzierte Modelle des übergebenen vollen Modells und gibt eine Vergleichstabelle mit den AICs zurück.
Königslachse
Prüfen der Annahmen am vollen Modell 1
par(mfrow = c(2,2)) plot(full_mod)
Königslachse
Prüfen der Annahmen am vollen Modell 2
Zusätzliche Plots: Residuen vs. jede X Variable
ChinookArg$residuals <- resid(full_mod) par(mfrow = c(1,2), mar = c(4,4,1,1)) boxplot(residuals ~loc, ChinookArg, main = "Residuals vs. loc"); abline(h = 0, lty = 2) plot(residuals ~ tl_log, ChinookArg, main = "Residuals vs. tl_log"); abline(h = 0, lty = 2)
Königslachse
Modellreduktion - Manuell
Wenn Interaktionsterme im Modell vorkommen, müssen auch die enthaltenen Variablen als Hauptterme vorkommen. Es darf also nicht tl_log oder loc im reduzierten Modell fehlen. Nur die Interaktion selbst kann aktuell rausgenommen werden:
Runde 1
red_mod <- update(full_mod, .~.-loc:tl_log) AIC(full_mod, red_mod)
df AIC full_mod 7 70.53732 red_mod 5 67.57480
- Das reduzierte Modell hat den niedrigeren AIC, es kann also der Interaktionsterm entfernt werden.
Runde 2
# Das reduzierte Modell wird das volle: full_mod <- red_mod # Einzelne X entfernen red_mod1 <- update(full_mod, .~.-loc) red_mod2 <- update(full_mod, .~.-tl_log) AIC(full_mod, red_mod1, red_mod2)
df AIC full_mod 5 67.57480 red_mod1 3 87.69126 red_mod2 4 224.06339
- Das volle Modell hat den besten AIC, daher wird nichts weiter entfernt.
Königslachse
Modellreduktion - automatisch
Built-in Funktion step()
step(full_mod, method = "backward")
Start: AIC=-249.3
w_log ~ loc + tl_log + loc:tl_log
Df Sum of Sq RSS AIC
- loc:tl_log 2 0.1011 10.965 -252.27
<none> 10.864 -249.31
Step: AIC=-252.27
w_log ~ loc + tl_log
Df Sum of Sq RSS AIC
<none> 10.965 -252.267
- loc 2 2.634 13.599 -232.151
- tl_log 1 34.175 45.140 -95.779
Call:
lm(formula = w_log ~ loc + tl_log, data = ChinookArg)
Coefficients:
(Intercept) locPetrohue locPuyehue tl_log
-8.1217 -0.2281 -0.4570 2.3049 
- Zum Schluss wird das ‘optimale’ Modell zusammengefasst.
- Auch hier wird das Modell mit beiden Variablen ohne Interaktion als ‘bestes’ Modell gewählt.
Königslachse
Validierung des finalen Modells 1
best_mod <- lm(w_log ~ loc + tl_log, ChinookArg) par(mfrow = c(2,2)) plot(best_mod)
Königslachse
Validierung des finalen Modells 2
Zusätzliche Plots: Residuen vs. jede X Variable
ChinookArg$residuals <- resid(best_mod) par(mfrow = c(1,2), mar = c(4,4,1,1)) boxplot(residuals ~loc, ChinookArg, main = "Residuals vs. loc"); abline(h = 0, lty = 2) plot(residuals ~ tl_log, ChinookArg, main = "Residuals vs. tl_log"); abline(h = 0, lty = 2)
Königslachse
Ergebnis 1 - Regressionkomponente
summary(best_mod)
Call:
lm(formula = w_log ~ loc + tl_log, data = ChinookArg)
Residuals:
Min 1Q Median 3Q Max
-0.60633 -0.19351 -0.00076 0.14352 1.61286
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -8.12169 0.56827 -14.292 < 2e-16 ***
locPetrohue -0.22813 0.08013 -2.847 0.00528 **
locPuyehue -0.45699 0.09231 -4.951 2.74e-06 ***
tl_log 2.30492 0.12563 18.347 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3186 on 108 degrees of freedom
Multiple R-squared: 0.8962, Adjusted R-squared: 0.8934
F-statistic: 311 on 3 and 108 DF, p-value: < 2.2e-16Königslachse
p-Werte für kategoriale Variablen
Problem der t-Tests und p-Werte
# Extraktion der Koeffizientenschätzung mit tidy() aus dem broom Paket broom::tidy(best_mod)
# A tibble: 4 x 5 term estimate std.error statistic p.value <chr> <dbl> <dbl> <dbl> <dbl> 1 (Intercept) -8.12 0.568 -14.3 1.19e-26 2 locPetrohue -0.228 0.0801 -2.85 5.28e- 3 3 locPuyehue -0.457 0.0923 -4.95 2.74e- 6 4 tl_log 2.30 0.126 18.3 5.72e-35
- Diese basieren auf dem Vergleich der einzelnen Faktorstufen mit der ‘Intercept’, also der Baseline-Faktorstufe.
- Möglicherweise unterscheiden sich diese nicht zur Baseline, aber untereinander!
Königslachse
Ergebnis 2 - ANOVA-Komponente
Um die ANOVA Tabelle zu erhalten, wo alle Faktorstufen miteinander verglichen werden, nutze die Funktion anova() für dein lm Objekt:
anova(best_mod)
Analysis of Variance Table
Response: w_log
Df Sum Sq Mean Sq F value Pr(>F)
loc 2 60.542 30.271 298.16 < 2.2e-16 ***
tl_log 1 34.175 34.175 336.61 < 2.2e-16 ***
Residuals 108 10.965 0.102
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1→ Dieser p-Wert für loc sollte im Bericht stehen, nicht die aus der summary Funktion!
Königlachse
Gleichung
coef(best_mod)
(Intercept) locPetrohue locPuyehue tl_log -8.1216871 -0.2281264 -0.4569901 2.3049203
Argentina: \(w.log_{A} = -8.12+2.30*tl.log\)
Petrohuye: \(w.log_{Pe}= (-8.12+-0.23)+2.30*tl.log = -8.35+2.30*tl.log\)
Puyehuye: \(w.log_{Pu} = (-8.12+-0.46)+2.30*tl.log = -8.58+2.30*tl.log\)
Königslachse
Visualisierung
ChinookArg %>% data_grid(
tl_log = seq_range(tl_log, 20),
loc = factor(levels(loc))
) %>%
add_predictions(best_mod) %>%
ggplot(aes(x = tl_log, group = loc,
colour = loc)) +
geom_line(aes(y = pred)) +
geom_point(data = ChinookArg,
mapping = aes(y = w_log)) +
ylab("w_log") +
scale_colour_manual(values =
c("forestgreen", "orange1",
"deeppink4")) +
theme_classic()
Ablaufprotokoll
Ablaufprotokoll
Beispiel für ‘backward selection’
-
Datenexploration:
- Ausreißer in einer der X Variablen? Variable entfernen oder transformieren.
- Ausreißer in Y? Behalten und hoffen, dass diese nicht zu einflussreich sind.
- Y nicht normalverteilt? Beibehalten und hoffen, dass die Residuen normalverteilt sind.
- Zusammenhang nicht linear? Variablen transformieren (entweder X, oder Y oder beide).
- X kollinear? Variable mit höchstem VIF-Wert entfernen oder nach biologischem Wissen entscheiden.
- Interaktionen? Evtl. explizit mit modellieren.
- Starte mit dem vollen Modell und beziehe Interaktionen nur mit ein, wenn Du vorher weißt, dass sie eine Rolle spielen könnten → nicht auf numerischen Output schauen → erst Annahmen überprüfen: wenn starke Heterogenität in den Residuen ggf. (weitere) Transformation.
- Führe die Modellauswahl durch, d.h. entferne schrittweise einzelne Variablen bis der AIC minimal ist.
- Annahmen des “optimalen” Modells prüfen.
- Wenn alles ok, weiter zu 8.
-
Wenn nicht:
- Residuen nicht normal? Keine Sorge, Homogenität ist wichtiger.
- Irgendwelche einflussreichen Beobachtungen? Variable(n) entfernen, Beobachtungen (ganze Zeilen) entfernen oder transformieren.
- Heterogenität, d.h. Muster im Streudiagramm Residuen vs. Fit?
- Abhängigkeit, d. h. Muster im Streudiagramm Residuen vs. jedes X (auch die im optimalen Modell ausgeschlossenen)?
- Passe das Modell erneut an und überprüfe wieder die Annahmen.
- Wenn diese nun erfüllt sind, betrachte das numerische Ergebnis: wie viel Variabilität wird erklärt? Ist der Einfluss von X signifikant? Wie stark ist der Einfluss?
- Stelle das Ergebnis in einer Tabelle und in einem Streudiagramm dar (einschließlich der Regressionslinie).
Muster in den Residuen und mögliche Lösungen
- Keine Homogenität & kein klares Muster in den Residuen gegen jedes X:
- Transformiere die Antwortvariable.
- Füge Interaktionen hinzu.
- Füge weitere erklärende Variablen hinzu.
- Anwendung generalisierter linearer Modelle (GLM) mit z.B. Poisson-Verteilung wenn die Daten Zählungen sind.
- Keine Homogenität & klares Muster in den Residuen gegen jedes X:
- Füge Interaktionen hinzu.
- Wenn die erklärende Variable, bei der ein Muster zu erkennen ist, aus dem endgültigen Modell ausgeschlossen wurde, wieder hinzufügen.
- Füge nicht-lineare Terme der erklärenden Variable hinzu (z. B. quadratische Terme).
- Anwendung generalisierter additiver Modelle (GAM).
- Homogenität gegeben, aber klares Muster in den Residuen gegen jedes X:
- Transformiere die entsprechenden erklärenden Variablen
- Wenn die erklärende Variable, bei der ein Muster zu erkennen ist, aus dem endgültigen Modell ausgeschlossen wurde, wieder hinzufügen.
- GAM anwenden.
Fragen?
Bei weiteren Fragen: saskia.otto(at)uni-hamburg.de
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License except for the borrowed and mentioned with proper source: statements.
Image on title and end slide: Section of an infrared satellite image showing the Larsen C ice shelf on the Antarctic Peninsula - USGS/NASA Landsat: A Crack of Light in the Polar Dark, Landsat 8 - TIRS, June 17, 2017 (under CC0 license)